home *** CD-ROM | disk | FTP | other *** search
/ Dr. Windows 3 / dr win3.zip / dr win3 / PROGRAMR / GSRC208A.ZIP / LSODA.C < prev    next >
C/C++ Source or Header  |  1993-08-17  |  62KB  |  2,337 lines
  1. /*
  2. From tam@dragonfly.wri.com Wed Apr 24 01:35:52 1991
  3. Return-Path: <tam>
  4. Date: Wed, 24 Apr 91 03:35:24 CDT
  5. From: tam@dragonfly.wri.com
  6. To: whitbeck@wheeler.wrc.unr.edu
  7. Subject: lsoda.c
  8. Cc: augenbau@sparc0.brc.uconn.edu
  9.  
  10.  
  11. I'm told by Steve Nichols at Georgia Tech that you are interested in
  12. a stiff integrator.  Here's a translation of the fortran code LSODA.
  13.  
  14. Please note
  15. that there is no comment.  The interface is the same as the FORTRAN
  16. code and I believe the documentation in LSODA will suffice.
  17. As usual, a free software comes with no guarantee.
  18.  
  19. Hon Wah Tam
  20. Wolfram Research, Inc.
  21. tam@wri.com
  22. */
  23.  
  24. /*
  25. I have done some additions to lsoda.c . These were mainly to fill the
  26. gap of some features that were available in the fortran code and were
  27. missing in the C version.
  28.  
  29. Changes are:
  30.          * all messages printed by lsoda routines will start with
  31.            a '\n' (instead of ending with one).
  32.          * xsetf() was added so that user can control printing
  33.            of messages by lsoda. xsetf should be called before
  34.            calling lsoda: xsetf(0) switches printing off
  35.                   xsetf(1) swithces printing on (default)
  36.            this implies one new global variable prfl (print flag).
  37.            xsetf(0) will stop *any* printing by lsoda functions.
  38.            Turning printing off means losing valuable information
  39.            but will not scramble your stderr output ;-)
  40.            This function is part of the original FORTRAN version.
  41.          * xsetf() and intdy() are now callable from outside the
  42.            library as assumed in the FORTRAN version.
  43.          * created lsoda.h that should be included in blocks
  44.            calling functions in lsoda's library.
  45.          * iwork5 can now have an extra value:
  46.               0 - no extra printing (default)
  47.               1 - print data on each method switch
  48.           ->  2 - print useful information at *each* stoda
  49.               step (one lsoda call has performs many stoda
  50.               calls) and also data on each method switch
  51.            note that a xsetf(0) call will prevent any printing even
  52.            if iwork5 > 0.
  53.          * hu, tn were made available as extern variables.
  54.          * the variable _ml has been renamed to _ml
  55.  
  56.                     Pedro Mendes
  57.  
  58.               eMail:      INTERNET: prm@aber.ac.uk
  59.               sMail:                Pedro Mendes,
  60.                         Dept. of Biological Sciences,
  61.                         University of Wales, Aberystwyth
  62.                         Dyfed, SY23 3DA,
  63.                         United Kingdom.
  64.  
  65.               phone: international: + 44 970 622353
  66.                       U.K.:     0970 622353
  67. */
  68.  
  69. #include <stdio.h>
  70. #include <math.h>
  71. /*#include <conio.h>
  72. #include "globals.h"
  73. #include "globvar.h"
  74. #include "datab.h"*/
  75. #ifndef _ZTC
  76. #include <malloc.h>
  77. #define mem_malloc malloc
  78. #define mem_free free
  79. #define mem_realloc realloc
  80. #else
  81. #include <stdlib.h>
  82. #define MEM_DEBUG 1
  83. #include "mem.h"
  84. #endif
  85.  
  86. #define max( a , b )  ( (a) > (b) ? (a) : (b) )
  87. #define min( a , b )  ( (a) < (b) ? (a) : (b) )
  88.  
  89. #define ETA 2.2204460492503131e-16
  90.  
  91. /* functions defined in other blocks */
  92. extern void daxpy(), dgefa(), dgesl(), dscal();
  93. extern double ddot();
  94. extern int idamax();
  95.  
  96. /* functions accessible from other blocks */
  97. void xsetf( int  mflag );
  98. void intdy();
  99.  
  100. /* internal functions */
  101. static void terminate( int *istate );
  102. static void terminate2( double *y, double *t );
  103. void freevectors( void );
  104. static void successreturn( double *y, double *t, int itask,
  105.                int ihit, double tcrit, int *istate );
  106. static void stoda( int neq, double *y, void (*f)() );
  107. static void endstoda( void );
  108. static void prja( int neq, double *y, void (*f)() );
  109. static void correction( int neq, double *y, void (*f)(), int *corflag,
  110.             double pnorm, double *del, double *delp,
  111.             double *told, int *ncf, double *rh, int *m );
  112. static void resetcoeff( void );
  113. static void methodswitch( double dsm, double pnorm, double *pdh, double *rh );
  114.  
  115. static void
  116.    solsy(),
  117.    cfode(),
  118.    ewset(),
  119.    scaleh(),
  120.    orderswitch(),
  121.    corfailure();
  122.  
  123. static double
  124.    vmnorm(), bnorm(), fnorm();
  125.  
  126. /* new static variable to enable optional printing (see comments) */
  127. static int prfl = 1;
  128.  
  129. /* newly added static variables */
  130. /* _ml was originally ml - changed to avoid conflict with one external ml */
  131. static int _ml, mu;
  132. static mord[3] = { 0, 12, 5 };
  133. static double sqrteta, *yp1, *yp2;
  134. static double sm1[13] = { 0., 0.5, 0.575, 0.55, 0.45, 0.35, 0.25,
  135.                   0.2, 0.15, 0.1, 0.075, 0.05, 0.025 };
  136.  
  137. /* variables accessible from other routines */
  138. double hu, tn, h, tsw, tolsf;
  139. int nst, nfe, nje, nqu, nq, imxer, mused, meth;
  140.  
  141. /* static variables for lsoda() */
  142. static double ccmax, el0, hmin, hmxi, rc;
  143. static int illin = 0, init = 0, mxstep, mxhnil, nhnil, ntrep = 0,
  144.        nslast, nyh, ierpj, iersl, jcur, jstart, kflag, l,
  145.        miter, maxord, maxcor, msbp, mxncf, n;
  146. static double pdnorm;
  147. static int ixpr = 0, jtyp, mxordn, mxords;
  148.  
  149. /* no static variable for prja(), solsy() */
  150. /* static variables for stoda() */
  151.  
  152. static double conit, crate, el[14], elco[13][14], hold, rmax,
  153.        tesco[13][4];
  154. static int ialth, ipup, lmax, meo, nslp;
  155. static double pdest, pdlast, ratio, cm1[13], cm2[6];
  156. static int icount, irflag;
  157.  
  158. /* static variable for block data */
  159. static int mesflg = 1;
  160.  
  161. /* static variables for various vectors and the Jacobian. */
  162. static double **yh, **wm, *ewt, *savf, *acor;
  163. static int *ipvt;
  164.  
  165. /*
  166.    The following are useful statistics.
  167.  
  168.    hu,
  169.    h,
  170.    tn,
  171.    tolsf,
  172.    tsw,
  173.    nst,
  174.    nfe,
  175.    nje,
  176.    nqu,
  177.    nq,
  178.    imxer,
  179.    mused,
  180.    meth
  181. */
  182.  
  183. /*
  184.    free allocated vectors
  185. */
  186.  
  187. void freevectors( void )
  188. {
  189.  int i, lenyh;
  190.  
  191.    lenyh = 1 + max( mxordn, mxords );
  192.    for( i=1; i<=lenyh; i++ ) mem_free( (void *) yh[i] );
  193.    for( i=1; i<=nyh; i++ ) mem_free( (void *) wm[i] );
  194.    mem_free( (void *) yh );
  195.    mem_free( (void *) wm );
  196.    mem_free( (void *) ewt );
  197.    mem_free( (void *) savf );
  198.    mem_free( (void *) acor );
  199.    mem_free( (void *) ipvt );
  200. }     /*   end freevectors   */
  201.  
  202. /*
  203.    Terminate lsoda due to illegal input.
  204. */
  205.  
  206. static void terminate( int *istate )
  207. {
  208.  if ( ( illin == 5 ) && prfl )
  209.  {
  210.   printf( "\nlsoda -- repeated occurrence of illegal input\n" );
  211.   printf(   "         run aborted.. apparent infinite loop" );
  212.  }
  213.  else
  214.  {
  215.   illin++;
  216.   *istate = -3;
  217.  }
  218. }         /*   end terminate   */
  219.  
  220.  
  221. /*
  222.    Terminate lsoda due to various error conditions.
  223. */
  224.  
  225.  
  226. static void terminate2( double *y, double *t )
  227. {
  228.  int i;
  229.  
  230.  yp1 = yh[1];
  231.  for ( i = 1 ; i <= n ; i++ ) y[i] = yp1[i];
  232.  *t = tn;
  233.  illin = 0;
  234.  freevectors();
  235.  return;
  236. }         /*   end terminate2   */
  237.  
  238.  
  239. /*
  240.    The following block handles all successful returns from lsoda.
  241.    If itask != 1, y is loaded from yh and t is set accordingly.
  242.    *Istate is set to 2, the illegal input counter is zeroed, and the
  243.    optional outputs are loaded into the work arrays before returning.
  244. */
  245.  
  246. static void successreturn( double *y, double *t, int itask,
  247.                int ihit, double tcrit, int *istate )
  248. {
  249.  int i;
  250.  
  251.  yp1 = yh[1];
  252.  for ( i = 1 ; i <= n ; i++ ) y[i] = yp1[i];
  253.  *t = tn;
  254.  if ( itask == 4 || itask == 5 )
  255.   if ( ihit )
  256.    *t = tcrit;
  257.  *istate = 2;
  258.  illin = 0;
  259.  freevectors();
  260. }   /*   end successreturn   */
  261.  
  262.  
  263. /*
  264.    In this version all memory allocated using malloc is freed upon exit.
  265.    Therefore *istate = 2 and *istate = 3 will not work.
  266. */
  267.  
  268. void lsoda( f, neq, y, t, tout, itol, rtol, atol, itask, istate,
  269.         iopt, jt, iwork1, iwork2, iwork5, iwork6, iwork7, iwork8,
  270.         iwork9, rwork1, rwork5, rwork6, rwork7 )
  271.  
  272. void (*f)();
  273. int neq, itol, itask, *istate, iopt, jt;
  274. int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
  275. double *y, *t, tout, *rtol, *atol;
  276. double rwork1, rwork5, rwork6, rwork7;
  277.  
  278. /*
  279.    If the user does not supply any of these values, the calling program
  280.    should initialize those untouched working variables to zero.
  281.  
  282.    _ml = iwork1
  283.    mu = iwork2
  284.    ixpr = iwork5
  285.    mxstep = iwork6
  286.    mxhnil = iwork7
  287.    mxordn = iwork8
  288.    mxords = iwork9
  289.  
  290.    tcrit = rwork1
  291.    h0 = rwork5
  292.    hmax = rwork6
  293.    hmin = rwork7
  294. */
  295.  
  296.  
  297. {
  298.    int mxstp0 = 500, mxhnl0 = 10;
  299.  
  300.    int i, i1, i2, iflag, kgo, lf0, lenyh, ihit;
  301.